%% This script uses the CHAMP EEF and ACE IEF to find the transfer function between them
clear;

%% JULIA data 
% load /home/mnair/projects/longp/alldays JULI_SEG ACE_SEG;
% JULI_SEG = JULI_SEG.*24.366*1e-3; %mV/m
% [Txy_short,F_short] = tfestimate(ACE_SEG,JULI_SEG,hanning(72),0,72,1/(5*60)); %Txy is agains calculated just to get
% [Cxy_short,F_short] = mscohere(ACE_SEG,JULI_SEG,hanning(72),0,72,1/(5*60));
% N_data = length(ACE_SEG) / 72;
% Err_short = sqrt( 1/(2*(N_data-1)) .* ( (1-Cxy_short)./Cxy_short ) ) .* abs(Txy_short);
% Err_short_lg = 0.434 * Err_short ./ abs (Txy_short);
 load /home/mnair/projects/longp/eef_data_mod.mat eef;

% I tried to see weather the TF magnitude plot can be reproduced by just
% two random variables. 
% std_eef = std(eef(:,6));
% mean_eef = mean(eef(:,6));
% eef_rand = mean_eef + randn(length(eef),1)*std_eef;
% eef(:,6) = eef_rand;
% eef(:,15) = eef_rand;
% 


load /home/mnair/projects/longp/OMNI_ELEC_new ace_all;

% Replace ace data with random signals

% L = isnan(ace_all(:,2));
% ace_all(L,2) = 0;
% std_ace = std(ace_all(:,2));
% mean_ace = mean(ace_all(:,2));
% ace_rand = mean_ace + randn(length(ace_all),1)*std_ace;
% ace_all(:,2) = ace_rand;
% ace_all(L,2) = 0 ;
% 

load /home/mnair/projects/longp/aplist.mat; 
fday_ap = fday_ap - datenum(2000,1,1);

% Initializations

ace_fday = floor(ace_all(:,1));
ap_lower_limit = 20; %lower limt if Ap
N_seg = 1;%the increasing counter for array JULI_SEG & ACE_SEG
N_data = 0;
phase_delay = 17;%minutes
mjd_date = datenum(2000,1,1);
plc = 'g';%plot color
tol = 1.8/24; % max time interval tolerated
min_time_length = 32/24; % minimum time length required in decimal days
np = 0;
lt_start = 10; % Ideal 9-15 (data permits 7-17)
lt_end = 14;
nd = 1;
coh_plot_lim = 0.1;
ace_interp_method = 'spline';
minimum_period_plot = 5;

des_int = 2; %desired sampling interval in hours
             % the sampling interval is set ~30 minutes above the
             % champ EQ revisit time (92 minutes)
len = floor( min_time_length *24 / des_int ); % This number is the  data length

ace_gap_tol = ( min_time_length + 6/24 ) * 0.1 ;
ace_gap_len_tol = floor ( ace_gap_tol * 1440/5 );
% the gap tolerence is 10% of the ace data length

% finding gaps in ACE data with gaps > 12 minutes
 
 nd = 1;
 np = 0;
 
 
 for i = 1: length(ace_all) - 1,

          
        if  isnan(ace_all(i,2)) && isnan(ace_all(i + 1 ,2 )) 
            
            np = np + 1;
        
        else
            if ( ace_all(i,1) - ace_all(i - np ,1) > 12/1440 )
                data_index_ace(nd,2) = i;
                data_index_ace(nd,1) = i - np;
                nd = nd + 1;
                               
            end;
                np = 0;
            
        end;
end;
    



% Calculating the effective electric field (LUHR EPS)
% E = 8IEF sqrt (64 + IEF^2). This has the effect of restructing amplitudes
% > 7 to within 8

L = ace_all(:,2) < 0;
temp = abs(ace_all(:,2));
Eeff = 8 * temp ./ sqrt (64 + temp.^2);
ace_all(:,2) = Eeff;
ace_all(L,2) = ace_all(L,2) * -1;



% down sampling and resampling to original interval

%The following script is to resample the ace data at a lower rate.
% The matlab resample apply a low-pass filter before resampling
% this allows for some remedy to the aliasing issues. This is a very
% important step. With out this sever aliasing issues can crrep into the
% spectra analysis.

L = isnan(ace_all(:,2));
y = interp1(ace_all(~L,1),ace_all(~L,2),ace_all(:,1),ace_interp_method);
ace_down = resample(y,1,des_int*60/5); % Checked the resampling time axis OK
ace_inter = interp1(ace_all(1:des_int*60/5:end,1),ace_down,ace_all(:,1));
ace_all(:,2) = ace_inter;


% make sure that the ACE data gaps > 12 min are filled with zeros using the
% data index determined above


for i = 1:length(data_index_ace),
    
    ace_all(data_index_ace(i,1):data_index_ace(i,2),2) = 0;
end;
    

 ace_gap_tol = ( min_time_length + 6/24 ) * 0.1 ;
 len = floor( min_time_length *24 / des_int ); % This number is the time length (also data length) in hours
 lt_incr = 1;          % An index used for counting the LT window numbers


 % iterate through LT
 
for lt_start = 9:0.5:16,
    
lt_end = lt_start + 0.5;
 nd = 1;
 np = 0;
% Divide CHAMP EEF data based on data gaps, LT etc

for i = 1: 36665,
    
    if (eef(i+1,1) - eef(i,1) <= tol && eef(i,4) > lt_start && eef(i,4) < lt_end ...
            && eef(i,1) - eef(i - np, 1) <= min_time_length )
        np = np + 1;
    else
       
        if ( eef(i,1) - eef(i - np ,1) >= min_time_length )
            data_index(nd,2) = i;
            data_index(nd,1) = i - np;
            np = 0;
            nd = nd + 1;
        else
            np=0;
            
        end;
    end;
    
end;

%for ap_lower_limit = 0:5:40,

                
% iterations
CHAMP_SEG = [];
ACE_SEG = [];
TIME_SEG = [];
N_data = 0;
N_seg = 1;
np = 0;
nd = 1;

for i = 1: size(data_index,1),
    
    
    L = ace_all(:,1) + phase_delay/(60*24) >= eef(data_index(i,1),1) - 3/24 ...
    & ace_all(:,1) + phase_delay/(60*24) <= eef(data_index(i,2),1)+ 3/24 ;
    
    if sum(L) > 0,
        ace_time = ace_all(L,1) +  phase_delay/(60*24) ;
        ace_data = ace_all(L,2);
        L = isnan(ace_data);
        if any(ace_data),
            if  max(diff([ace_time(1) ; ace_time(~L) ; ace_time(end)] )) < ace_gap_tol && ...
                    sum(L) < ace_gap_len_tol,
                
                if abs(mean(diff(ace_time))-0.0035) <= 1e-004, %Use this with ACE 5 min averages
                    %if abs(mean(diff(ace_time))-2/24) <= 1e-004, %Use this with ACE 2 hours averages
                    dummy = eef(data_index(i,1):data_index(i,2), 6 ) - eef(data_index(i,1):data_index(i,2), 15 )/1e3;
                    
                    % resampling the CHAMP data (This step is not changing
                    % the results much. )
%                     [n,d] = rat(median(diff(eef(data_index(i,1):data_index(i,2),1)))/(des_int/24));
%                     y = resample(dummy,  n,  d);  % Now resample it
%                     t2 = (0:(length(y)-1)) * d / (n * 1/ median(diff(eef(data_index(i,1):data_index(i,2),1))) );
                    
                    time_axis_desired = eef(data_index(i,1), 1 ):des_int/24:eef(data_index(i,2), 1 );
                    CHAMP_EEF = interp1(eef(data_index(i,1):data_index(i,2),1), dummy,time_axis_desired);
                    %CHAMP_EEF = interp1(t2 + eef(data_index(i,1),1), y, time_axis_desired);
                    
                    ACED = interp1(ace_time, ace_data, time_axis_desired);
                    
                    L = isnan(ACED);
                    ACED(L) = [];
                    
                    L = fday_ap >= eef(data_index(i,1), 1 )& fday_ap <= eef(data_index(i,2), 1 );
                    mean_ap = mean(ap(L));
                    
                    if sum(isnan(ACED)) < 1 & mean_ap >= ap_lower_limit & length(ACED) >= len & length(CHAMP_EEF) >= len;% & mean_ap <=20,
                        
                        
                        CHAMP_SEG(N_seg:N_seg+(len-1)) = CHAMP_EEF(1:len);
                        ACE_SEG(N_seg:N_seg+(len-1)) = ACED(1:len);
                        TIME_SEG(N_seg:N_seg+(len-1)) = time_axis_desired(1:len);
                        % IMF_BZ_SEG(N_seg:N_seg+71) = IMFBZ(1:72);
                        N_seg = N_seg+len;
                        N_data = N_data+1;
                        
                        
                    end;
                else,
                    %fprintf('Day %d has some missing time stamp\n', Julia_W(i).fday);
                end;
            end;
        end;
    end;
end;

fprintf('Done ! N_data = %3d , Length = %d hours, LT = %d-%d, AP limit = %d color = %s, interp = %s\n', ...
    N_data, min_time_length*24, lt_start,lt_end, ap_lower_limit, plc, ace_interp_method);

% Coherence, phase and tranfer function
%JULI_SEG = JULI_SEG.*24.366*1e-3; %mV/m
figure(1);
[Cxy_long,F] = mscohere(ACE_SEG, CHAMP_SEG*1e3,hanning(len),0,len,1/(des_int*3600)); %1/(5*60) = sampling frequency in Hz )
[Pxy,F] = cpsd(CHAMP_SEG*1e3,ACE_SEG,hanning(len),0,len,1/(des_int*3600)); %
[Pxx,F] = pwelch(CHAMP_SEG*1e3,hanning(len),0,len,1/(des_int*3600));
[Pyy,F] = pwelch(ACE_SEG,hanning(len),0,len,1/(des_int*3600));
phase = angle(Pxy);
[Txy_long,F_long] = tfestimate(ACE_SEG,CHAMP_SEG*1e3, hanning(len),0,len,1/(des_int*3600));
%tf = conj(Txy');
L = Cxy_long > coh_plot_lim & (1./F_long)./(3600) >= minimum_period_plot ;

Err_long = sqrt( 1/(2*(N_data-1)) .* ( (1-Cxy_long)./Cxy_long ) ) .* abs(Txy_long);

Err_long_lg = 0.434 * Err_long ./ abs (Txy_long);

% errorbar(log10((1./F_long(L))./(3600)), log10(abs(Txy_long(L))),(Err_long_lg(L)),(Err_long_lg(L)),plc, 'linewidth',2);
% hold on;
% axis([ -inf   inf   -3   -1]);

for i = 2:length(Cxy_long),
    fprintf('%06.3f ',(1./F_long(i))./(3600));
end;
fprintf('\n');
for i = 2:length(Cxy_long),
    fprintf('%06.3f ', Cxy_long(i)); 
end;
fprintf('\n');
%%
%errorbar(log10((1./F_short(2:end))./(3600)), log10(abs(Txy_short(2:end))),log10(Err_short(2:end))/2,log10(Err_short(2:end))/2,'g');

%errorbar(log10((1./F_short(2:end))./(3600)), log10(abs(Txy_short(2:end))),Err_short_lg(2:end)/2,Err_short_lg(2:end)/2,'k');
%errorbar(log10((1./F_short(2:end))./(3600)), log10(abs(Txy_short(2:end))),(Err_short_lg(2:end)),(Err_short_lg(2:end)),'b');

title_str = sprintf(' Proc%3d-%dLT%d-%dAP%d%s\n',N_data, min_time_length*24, lt_start,lt_end, ap_lower_limit, ace_interp_method);
%title(title_str);

% saveas(gcf, ['/home/mnair/projects/longp/pic/' title_str],'png');
% 
% close all;
% 

lt_axis(lt_incr,:) = (lt_start + lt_end ) / 2;
cohdata(lt_incr,:) = Cxy_long;
pxxdata(lt_incr,:) = Pxx;
pyydata(lt_incr,:) = Pyy;
pxydata(lt_incr,:) = Pxy;
lt_incr = lt_incr + 1;
end;
